Statistische Modellierung mit R

Ein Workshop für die PH Zürich

Veröffentlichungsdatum

17. Oktober 2025

Methodology: The Science before Statistics

Jede statistische Modellierung gewinnt an Aussagekraft, je umfassender sie die inhaltliche Fragestellung abzubilden im Stande ist. Um aus der riesigen Fülle an Optionen geeignet und zielgerichtet auswählen zu können sind die folgenden Unterscheidungen oft sehr hilfreich.

Erkenntnisinteressen

Ganz grundlegend kann a priori das Erkenntnisinteresse von Studien in die folgenden vier Kategorien unterschieden werden:

Erkenntisinteressen nach (Döring und Bortz 2016).
Deskriptiv Explorativ Explanativ Prädiktiv
populationsbeschreibend hypothesengenerierend hypothesenprüfend Datenpunkte vorhersagend oder imputierend
Bei welchem Anteil 15-Jähriger in Deutschland handelt es sich um funktionale Analphabet:innen? Was sind potentielle Ursachen für genderbezogene Disparitäten im Analphabetismus? Sind 15-jährige Jungen häufiger Analphabeten als 15-jährige Mädchen? Mit welchen Variablen können Schüler:innen at risk erfolgreich identifiziert werden?

Gütekriterien wiss. Erkenntnis nach Campbell (1957)

Für ein erfolgreiches Studiendesign und die anschließende statistische Analyse ist es sehr wertvoll sich vorab über Schwerpunkte besonders gewünschter Aspekte wissenschaftlicher Güte Gedanken zu machen. Insbesondere über die Unterkriterien Methodischer Strenge:

  • Konstruktvalidität (Inwiefern ist die Interpretation der Messwerte angemessen?)
  • Interne Validität (Inwiefern sind Assoziationen von unabhängiger [beeinflussender] und abhängiger [beeinflusster] Variabler als kausale Effekte interpretierbar?)
  • Externe Validität (Inwiefern können die Schlussfolgerungen der Studie verallgemeinert werden?)
  • Statistische Validität (Wie robust und angemessen sind die verwendeten statistischen Verfahren?)

Deskriptiv- und Inferenzstatistik

Deskriptive Statistik beschreibt vorliegende Daten (z.B. mit Effektstärken), während Inferenzstatistik Aussagen über den die Daten generierenden Mechanismus trifft. Beide können »eher einfach« oder »hoch komplex« sein und oftmals stehen sie in einem synergetischen Verhältnis (siehe Abbildung 1).

Abbildung 1: Synergetisches Verhältnis von deskriptiver Statistik und Inferenzstatistik

Bayesianisches und Frequentistisches Schätzen und Testen

Eine sehr heuristische Klassifikation inferenzstatistischer Verfahren stellt die Unterscheidung von statistischer Schätzung und Testung dar:

(Inferenzstatistische) Schätzungen (estimation with quantified uncertainty) treffen anhand von Stichproben Aussagen über Parameter der Grundgesamtheit (Population) aus der die Stichprobe gezogen wurde.1

1 Bspw.: »Mit 96%iger Wahrscheinlichkeit liegt die Analphabetismusinzidienz von 15-Jährigen in Deutschland zwischen .08 und .12«

(Inferenzstatistische) Hypothesentests bewerten anhand von Stichprobendaten die Gültigkeit von Hypothesen in der Grundgesamtheit (Population) aus der die Stichprobe gezogen wurde.2

2 Bspw.: »Nimmt man an, dass sich die Analphabetismusinzidienz von 15-Jährigen in Deutschland zwischen 2021 und 2025 nicht geändert hat, beträgt die Wahrscheinlichkeit der vorliegenden Daten p = .032«

Diese beiden Verfahren können sowohl im Rahmen der frequentistischen Statistik als auch der bayesianischen Statistik angewendet werden. Die folgende Tabelle gibt einen Überblick über die wichtigsten Werkzeuge:

Frequentistische Statistik Bayesianische Statistik
Parameterschätzung Konfidenzintervalle Posterior Distributions
Hypothesentest p-Werte Bayes Faktoren & ROPE-CrI Procedure

Hypothesenarten

Bayesianische wie frequentistischen Hypothesentests können unterschiedliche Arten von Hypothesen zugrunde gelegt werden:

  • Punkthypothesen setzen Parameter gleich einer reellen Zahl; etwa \(H_0\text{: } \delta = 0\)
  • Äquivalenzhypothesen nehmen Parameter in einem reellen Intervall an; etwa \(H_0\text{: } \delta \not\in\ [-.3, .3]\)
  • Informative Hypothesen nehmen eine Ordnungsrelation mehrerer Parameter an; etwa \(\mu_{\text{Baseline}} < \mu_{\text{Imaginary Pill}} < \mu_{\text{Blinded Placebo}}\) (Buergler u. a. 2023)

Die Art der (falsifizierten) Hypothese entscheidet wesentlich stärker über den Informationsgehalt eines Hypothesentests als die Entscheidung für das frequentistische oder bayesianische Paradigma (Hoijtink 2012).

Dies ist am leichtesten anhand der Nullhypothese nachvollziehbar. Wird etwa die Nullhypothese \(H_0\text{: } \delta = 0\) verworfen, wird entsprechend die Alternativhypothese \(H_A\text{: } \delta \neq 0\) angenommen. Diese enthält aber quasi keine Information, da sie nur mit einer einzigen Beobachtung (d = 0.000000 …) verworfen werden kann und im kritischen Rationalismus gilt, dass eine Aussage umso mehr Information enthält, umso leichter sie verworfen werden kann (Döring und Bortz 2016).

Äquivalenzhypothesen können sowohl frequentistisch (z.B. TOAST-Prozedur in R und JASP, Lakens 2017) wie bayesianisch (z.B. ROPE-Ansatz Kruschke 2015) getestet werden. Für das Testen informativer Hypothesen liegen bayesianische Methoden in (u.a.) JASP und R vor (z.B. {bain}, Gu u. a. 2019) sowie in den frequentistischen R-Paketen {restriktor} (Vanbrabant 2020) und {ic.infer} (Grömping 2010).

Grundlagen der Regressionsanalyse

Regressionsanalysen sind ein sehr mächtiges Werkzeug um Zusammenhänge oder Unterschiede zwischen/in Variablen zu modellieren. Innerhalb der Regressionsanalyse kann sowohl geschätzt als auch getestet werden und dies jeweils bayesianisch oder frequentistisch.

Die Grundidee der lineare Regression ist in diesem interaktiven Applet veranschaulicht. Eine Abhängige Variable \(y_i\) wird also in der Regressionsanlyse als normalverteilt mit bedingtem Erwartungswert \(b_0 + b_1*x_i\) und \(\sigma\) angenommen:

\[y_i \sim \mathcal{N}\left( b_0 + b_1*x_i, \sigma \right)\]

Datengrundlage

Das wollen wir in den berühmten Daten aus dem STAR-Experiment (Student/Teacher Achievement Ratio) illustrieren. In diesem Experiment wurden Schüler:innen und Lehrer:innen zufällig auf Klassen mit kleiner (13-17 Schüler:innen pro Lehrer:in) und großer (22-25 Schüler:innen pro Lehrer:in) Klassengröße zugeteilt. Zusätzlich gab es Klassen großer Größe, die von einer:m ausgebildeten Hilfslehrer:in unterstützt wurde. Die Schüler:innen wurden über mehrere Jahre getestet.

Wir können die Daten herunterladen oder direkt programmatisch importieren.

library(readr) # für den data import
library(dplyr)    # für das datawrangling
data_star <-
  readr::read_csv(
    paste0(
      "https://raw.githubusercontent.com/sammerk/",
      "Zuerich_Statistics/master/data/star.csv"
    )
  )

Der Satzsatz sieht wie folgt aus:

glimpse(data_star)
Rows: 26,796
Columns: 18
$ id      <dbl> 100017, 100028, 100045, 100045, 100045, 100064, 100070, 100070…
$ sch     <dbl> 28, 52, 41, 41, 41, 64, 40, 40, 22, 22, 22, 17, 17, 3, 3, 3, 3…
$ gr      <chr> "K", "K", "1", "2", "3", "K", "K", "1", "K", "1", "2", "K", "1…
$ cltype  <chr> "small", "reg", "small", "small", "small", "small", "reg", "re…
$ hdeg    <chr> "BS/BA", "MS/MA/MEd", "BS/BA", "MS/MA/MEd", "BS/BA", "BS/BA", …
$ clad    <chr> "1", "1", "1", "APPR", "1", "1", "APPR", "1", "PROB", "1", "NO…
$ exp     <dbl> 3, 12, 20, 15, 5, 19, 2, 5, 9, 9, 2, 20, 16, 8, 23, 21, 14, 6,…
$ trace   <chr> "B", "W", "W", "B", "W", "W", "W", "W", "B", "B", "B", "W", "W…
$ read    <dbl> 476, 410, 507, 575, 610, 430, 495, 629, 418, 524, 627, 421, 45…
$ math    <dbl> 602, 444, 572, 620, 655, 484, 576, 592, 489, 567, 603, 459, 46…
$ ses     <chr> "F", "N", "N", "N", "N", "N", "F", "F", "F", "F", "F", "N", "N…
$ schtype <chr> "inner", "suburb", "suburb", "suburb", "suburb", "rural", "rur…
$ sx      <chr> "F", "F", "M", "M", "M", "F", "F", "F", "F", "F", "F", "M", "M…
$ eth     <chr> "B", "W", "W", "W", "W", "W", "W", "W", "B", "B", "B", "W", "W…
$ birthq  <chr> "1980:1", "1980:2", "1980:2", "1980:2", "1980:2", "1980:1", "1…
$ birthy  <dbl> 1980, 1980, 1980, 1980, 1980, 1980, 1979, 1979, 1980, 1980, 19…
$ yrs     <dbl> 0, 0, 1, 2, 3, 0, 0, 1, 0, 1, 2, 0, 1, 0, 2, 3, 0, 0, 0, 0, 1,…
$ tch     <dbl> 478, 893, 698, 701, 706, 1102, 679, 684, 351, 359, 365, 267, 2…

Die Variablen dieses Datensatzes sind die folgenden:

  • id: a factor - student id number
  • sch: a factor - school id number
  • gr: grade - an ordered factor with levels K < 1 < 2 < 3
  • cltype: class type - a factor with levels small, reg and reg+A. The last level indicates a regular class size with a teachers aide.
  • hdeg: highest degree obtained by the teacher - an ordered factor with levels ASSOC < BS/BA < MS/MA/MEd < MA+ < Ed.S < Ed.D/Ph.D
  • clad: career ladder position of the teacher - a factor with levels NOT APPR PROB PEND 1 2 3
  • exp: a numeric vector - the total number of years of experience of the teacher
  • trace: teacher’s race - a factor with levels W, B, A, H, I and O representing white, black, Asian, Hispanic, Indian (Native American) and other
  • read: the student’s total reading scaled score
  • math: the student’s total math scaled score
  • ses: socioeconomic status - a factor with levels F and N representing eligible for free lunches or not eligible
  • schtype: school type - a factor with levels inner, suburb, rural and urban
  • sx: student’s sex - a factor with levels M F
  • eth: student’s ethnicity - a factor with the same levels as trace
  • birthq: student’s birth quarter - an ordered factor with levels 1977:1 < … < 1982:2
  • birthy: student’s birth year - an ordered factor with levels 1977:1982
  • yrs: number of years of schooling for the student - a numeric version of the grade gr with Kindergarten represented as 0. This variable was generated from gr and does not allow for a student being retained.
  • tch: a factor - teacher id number

Um nicht gleich mit Mehrebenenregression starten zu müssen, erstellen wir als erstes einen Subdatensatz, der nur eine:n Schüler:in pro Lehrer:in enthält

data_star_sub <-
  data_star %>%
  group_by(tch) %>%
  sample_n(1) %>%
  ungroup()

»Effekte« der Klassenstufe

Einfache lineare Regression

Zunächst wollen wir »Effekte« der Klassenstufe modellieren. Ich empfehle dingend Daten immer erst zu visualisieren und dann statistisch zu modellieren.

library(ggplot2) # plots
Warning: package 'ggplot2' was built under R version 4.4.1
library(ggforce) # sina plots

ggplot(data_star_sub, aes(gr, math)) +
  geom_jitter() +
  theme_minimal()
Warning: Removed 120 rows containing missing values or values outside the scale range
(`geom_point()`).

Die unabhängige Variable ist hier noch nicht intervallskaliert. Dies können wir wie folgt ändern.

# Das Schuljahr K zu 0 rekodieren
data_star_sub <-
  data_star_sub %>%
  mutate(
    grade = case_when(
      gr == "K" ~ 0,
      gr == "1" ~ 1,
      gr == "2" ~ 2,
      gr == "3" ~ 3
    )
  )
# Anschließend nochmals visualisieren
ggplot(data_star_sub, aes(grade, math)) +
  geom_jitter() +
  theme_minimal()
Warning: Removed 120 rows containing missing values or values outside the scale range
(`geom_point()`).

Dann kann eine erste Regressionsgerade modelliert werden.

Frequentistische Punktschätzung

Mit der Funktion lm() werden mithilfe kleinster Quadrate \(b_0\) und \(b_1\) geschätzt.

mod0 <- lm(math ~ grade, data = data_star_sub)
mod0

Call:
lm(formula = math ~ grade, data = data_star_sub)

Coefficients:
(Intercept)        grade  
     485.16        46.64

Frequentistischer Nullhypothesensignifikanztest

Die Standardinferenzstatistik zu einem solchen Modell ist sind p-Werte für die Nullhypothesen \[H_0:\; b_0 = 0\] \[H_0:\; b_1 = 0\]

summary(mod0)

Call:
lm(formula = math ~ grade, data = data_star_sub)

Residuals:
     Min       1Q   Median       3Q      Max 
-197.162  -28.080   -1.162   25.920  142.559 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  485.162      1.980  245.04   <2e-16 ***
grade         46.639      1.083   43.08   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 42.55 on 1265 degrees of freedom
  (120 observations deleted due to missingness)
Multiple R-squared:  0.5946,    Adjusted R-squared:  0.5943 
F-statistic:  1856 on 1 and 1265 DF,  p-value: < 2.2e-16

Frequentistische Konfidenzintervalle

confint(mod0)
                2.5 %    97.5 %
(Intercept) 481.27755 489.04613
grade        44.51531  48.76346

Bayesianische Posterior Distributions

library(MCMCpack) #für bayesianische Schätzung
mod1 <- MCMCregress(math ~ grade, data = data_star_sub)

plot(mod1)

summary(mod1)

Iterations = 1001:11000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 10000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

               Mean     SD Naive SE Time-series SE
(Intercept)  485.17  1.980  0.01980        0.01980
grade         46.64  1.095  0.01095        0.01095
sigma2      1813.72 72.286  0.72286        0.72286

2. Quantiles for each variable:

              2.5%    25%     50%     75%   97.5%
(Intercept)  481.3  483.8  485.18  486.50  489.06
grade         44.5   45.9   46.62   47.36   48.78
sigma2      1677.0 1764.0 1811.87 1861.82 1959.01

Während die Densityplots den Posterior der Schätzung charakterisieren, erhält man mit summary() auch die 95%CrI-Intervalle.

Bayes Factor

Mit einem Bayes Factor testet man im Falle unseres Models wie wahrscheinlich die vorliegenden Daten sind geben ein Intercept-Only-Modell vs. unser Model mod1. Da die Funktion lmBF() im Paket {BayesFactor} allerdings Daten ohne Missings erwartet, müssen diese noch gefiltert werden.

library(BayesFactor)
Warning: package 'Matrix' was built under R version 4.4.1
mod2 <- lmBF(
  formula = math ~ grade,
  data = data_star_sub %>%
    filter(!is.na(math))
)
Warning: data coerced from tibble to data frame

Standardisierte Regression

Die in mod0 erhaltenen Regressionskoeffizienten konnten wir nutzen um zu interpretieren welche Punktzahl durchschnittlich im Kindergarten vorliegt (\(b_0\)) und was der durchschnittliche Anstieg in einem Jahr ist (\(b_1\)). Diese Zahlen sind insbesondere dann gut als (nicht-standardisierte) Effektstärke interpretierbar, wenn man mit der Skala vertraut ist. Eine Alternative dazu ist die z-Standardisierung der abhängigen Variable. Um dies zu veranschaulichen betrachten wir Effekte des SES (free lunch) auf die Mathematikleistung in der dritten Klasse. Die Rohdataten sehen wie folgt aus:

data_star_sub %>%
  filter(grade == 3 & !is.na(ses)) %>%
  ggplot(aes(ses, math)) +
  geom_violin() +
  geom_sina() +
  stat_summary(
    fun.data = mean_sdl,
    fun.args = list(mult = 1),
    geom = "pointrange",
    colour = "#1bbc9d"
  ) +
  theme_minimal()
Warning: Removed 32 rows containing non-finite outside the scale range
(`stat_ydensity()`).
Warning: Removed 32 rows containing non-finite outside the scale range
(`stat_sina()`).
Warning: Removed 32 rows containing non-finite outside the scale range
(`stat_summary()`).

z-standardisieren wir die abhängige Variable bleiben die Verteilungsformen gleich, lediglich die Skalierung ändert sich.

data_star_sub %>% 
  filter(grade == 3 & !is.na(ses)) %>% 
  ggplot(aes(ses, scale(math))) +
  geom_violin() +
  geom_sina() +
  stat_summary(
    fun.data = mean_sdl,
    fun.args = list(mult = 1),
    geom = "pointrange",
    colour = "#1bbc9d"
  ) +
  theme_minimal()
Warning: Removed 32 rows containing non-finite outside the scale range
(`stat_ydensity()`).
Warning: Removed 32 rows containing non-finite outside the scale range
(`stat_sina()`).
Warning: Removed 32 rows containing non-finite outside the scale range
(`stat_summary()`).

Der folgende Code kodiert automatisch die ses Variable als 0, 1-Variable um.

mod3 <-
  lm(
    scale(math) ~ ses,
    data = data_star_sub %>%
      filter(grade == 3 & !is.na(ses))
  )

Die Koeffizienten passend sehr gut zur Abbildung:

\[ \operatorname{\widehat{scale(math)}} = -0.27 + 0.49(\operatorname{ses}_{\operatorname{N}}) \]

-0.27 ist der Mittelwert der Gruppe mit Free Lunch, 0.49 die Differenz der beiden Gruppenmittelwerte gemessen in Standardabweichungen und damit wie ein Cohen’s \(d\) interpretiertbar.

library(effectsize)
Warning: package 'effectsize' was built under R version 4.4.1
cohens_d(
  math ~ ses,
  data = data_star_sub %>% filter(grade == 3 & !is.na(ses)),
  pooled_sd = F
)
Warning: Missing values detected. NAs dropped.
Cohen's d |         95% CI
--------------------------
-0.51     | [-0.74, -0.27]

- Estimated using un-pooled SD.

Multiple Regression

Die multiple Regression erweitert die Idee der einfachen lineare Regression auf mehrere unabhängige Variablen:

\[y_i \sim \mathcal{N}\left( b_0 + b_1*x_{1i} + b_2*x_{2i}, \sigma \right)\]

Informative Hypothesentests in der multiplen Regression

# Informative Hypothese
library(bain)

## ANOVA als Regression

data_star_sub <- data_star_sub %>% 
    mutate(classtype = ifelse(cltype == "reg+A",
                                                        "regA",  cltype))

mod5 <- lm(math ~ classtype - 1, data = data_star_sub) # Intercept mit `-1` entfernt
mod5

Call:
lm(formula = math ~ classtype - 1, data = data_star_sub)

Coefficients:
  classtypereg   classtyperegA  classtypesmall  
         544.3           555.4           558.3  
bain(mod5, "classtypereg = classtyperegA = classtypesmall; classtypereg < classtyperegA < classtypesmall")
Bayesian informative hypothesis testing for an object of class lm (continuous predictors):

   Fit   Com   BF.u  BF.c   PMPa  PMPb  PMPc 
H1 0.000 0.000 4.095 4.095  0.478 0.428 0.461
H2 0.734 0.164 4.477 14.064 0.522 0.468 0.504
Hu                                0.104      
Hc 0.266 0.836 0.318                    0.036

Hypotheses:
  H1: classtypereg=classtyperegA=classtypesmall
  H2: classtypereg<classtyperegA<classtypesmall

Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement. PMPa contains the posterior model probabilities of the hypotheses specified. PMPb adds Hu, the unconstrained hypothesis. PMPc adds Hc, the complement of the union of the hypotheses specified.

Schätzung der eines Cohen’s d zwischen großer und kleiner Klasse in math

data_star %>% 
    filter(gr == 3 & !is.na(cltype) & cltype != "reg+A") %>%
    cohens_d(
      math ~ cltype,
      data = .,
      pooled_sd = F
    )
Warning: Missing values detected. NAs dropped.
Cohen's d |         95% CI
--------------------------
-0.17     | [-0.24, -0.11]

- Estimated using un-pooled SD.

Literatur

Buergler, Sarah, Dilan Sezer, Niels Bagge, Irving Kirsch, Cosima Locher, Claudia Carvalho, und Jens Gaab. 2023. „Imaginary Pills and Open-Label Placebos Can Reduce Test Anxiety by Means of Placebo Mechanisms Scientific Reports“. Scientific Reports 13 (1): 2624. https://doi.org/10.1038/s41598-023-29624-7.
Campbell, Donald T. 1957. „Factors Relevant to the Validity of Experiments in Social Settings“. Psychological Bulletin 54 (4): 297–312. https://doi.org/10.1037/h0040950.
Döring, Nicola, und Jürgen Bortz. 2016. Forschungsmethoden und Evaluation in den Sozial- und Humanwissenschaften. 5., vollst. Berlin, Heidelberg: Springer.
Grömping, Ulrike. 2010. „Inference with Linear Equality and Inequality Constraints Using R: The Package Ic.infer“. Journal of Statistical Software 33 (Februar): 1–31. https://doi.org/10.18637/jss.v033.i10.
Gu, Xin, Herbert Hoijtink, Joris Mulder, und Yves Rosseel. 2019. „Bain: A Program for Bayesian Testing of Order Constrained Hypotheses in Structural Equation Models“. Journal of Statistical Computation and Simulation 89 (8): 1526–53. https://doi.org/10.1080/00949655.2019.1590574.
Hoijtink, Herbert. 2012. Informative Hypotheses: Theory and Practice for Behavioral and Social Scientists. Chapman & Hall/CRC Statistics in the Social and Behavioral Sciences Series. Boca Raton: CRC.
Kruschke, John K. 2015. Doing Bayesian Data Analysis: A Tutorial with R, JAGS, and Stan. 2nd Aufl. Academic Press.
Lakens, Daniël. 2017. „Equivalence Tests: A Practical Primer for t Tests, Correlations, and Meta-Analyses“. Social Psychological and Personality Science 8 (4): 355–62. https://doi.org/10.1177/1948550617697177.
Vanbrabant, L. 2020. Restriktor: Constrained Statistical Inference. https://cran.r-project.org/package=restriktor.